//-------------------------------------------------------------------
// Purpose: Compare the speed of two algoritms for finding the closest
//          pair of points in an array of points.  The brute force 
//          method compares every point to every other point.  The 
//          recursive method splits the points based on x coordinate,
//          finds the closest pair of points in each half, and then
//          combines the results to get the solution.
//
// Author:  John Gauch
//-------------------------------------------------------------------

#include <cstdlib>
#include <cmath>
#include <iostream>
using namespace std;

const int MAX_SIZE = 10000;
int sqrt_count = 0;

//-------------------------------------------------------------------
// Simple (x,y) point class.
//-------------------------------------------------------------------
class Point 
{
public:
   int x,y;
};

//-------------------------------------------------------------------
// Array of point class.
//-------------------------------------------------------------------
class PointArray
{
public:
   void random(int range);
   void print();
   void xsort();
   void ysort();
   Point point[MAX_SIZE];
   int size;
};

//-------------------------------------------------------------------
// Assign random values to array of points.
//-------------------------------------------------------------------
void PointArray::random(int count)
{
   size = count;
   int range = size * 10;
   for (int i=0; i<size; i++)
   {
      point[i].x = rand() % range;
      point[i].y = rand() % range;
   }
}

//-------------------------------------------------------------------
// Print array of points.
//-------------------------------------------------------------------
void PointArray::print()
{
   for (int i=0; i<size; i++)
      cout << i << " " << point[i].x << " " << point[i].y << endl;
}

//-------------------------------------------------------------------
// Sort array of points in x direction
//-------------------------------------------------------------------
void PointArray::xsort()
{
   for (int i=0; i<size; i++)
   for (int j=i+1; j<size; j++)
   if ((point[i].x > point[j].x) ||
      ((point[i].x == point[j].x) && (point[i].y > point[j].y)))
   {
      Point temp = point[i];
      point[i] = point[j];
      point[j] = temp;
   }
}

//-------------------------------------------------------------------
// Sort array of points in y direction
//-------------------------------------------------------------------
void PointArray::ysort()
{
   for (int i=0; i<size; i++)
   for (int j=i+1; j<size; j++)
   if ((point[i].y > point[j].y) ||
      ((point[i].y == point[j].y) && (point[i].x > point[j].x)))
   {
      Point temp = point[i];
      point[i] = point[j];
      point[j] = temp;
   }
}

//-------------------------------------------------------------------
// Brute force algorithm to find closest pair of points.
//-------------------------------------------------------------------
void closest_pair(PointArray &pts,
                  Point &p1, Point &p2, double &dist)
{
   dist = MAX_SIZE;
   for (int i=0; i<pts.size; i++)
   for (int j=i+1; j<pts.size; j++)
   {
      int dx = pts.point[i].x - pts.point[j].x;
      int dy = pts.point[i].y - pts.point[j].y;
      double temp = sqrt(dx*dx + dy*dy);
      sqrt_count++;
      if (temp < dist)
      {
         dist = temp;
         p1 = pts.point[i];
         p2 = pts.point[j];
      }
   }
}

//-------------------------------------------------------------------
// Split y array into two parts based on the midpoint value in x array.
//-------------------------------------------------------------------
void tricky_split(PointArray &xpts, int xstart, int xend, 
                  PointArray &ypts, int ystart, int yend, 
                  PointArray &left, PointArray &right, 
                  int &xmid, int &xvalue, int &yvalue)
{
   // Select midpoint of xpts array
   xmid = (xstart + xend) / 2;
   xvalue = xpts.point[xmid].x;
   yvalue = xpts.point[xmid].y;

   // Get points on left and right of xvalue
   left.size = 0;
   right.size = 0;
   for (int i = ystart; i < yend; i++)
   {
      if ((ypts.point[i].x < xvalue) ||
         ((ypts.point[i].x == xvalue) && (ypts.point[i].y < yvalue)))
         left.point[left.size++] = ypts.point[i];
      else 
         right.point[right.size++] = ypts.point[i];
   }
}

//-------------------------------------------------------------------
// Merge closest pair results by comparing points near midpoint line.
//-------------------------------------------------------------------
void tricky_merge(int xvalue, PointArray &pts, int start, int end, 
                  Point &p1, Point &p2, double &dist)
{
   // Get points on left of xvalue
   for (int i = start; i < end; i++)
   {
      int dx = pts.point[i].x - xvalue;
      if ((-dist <= dx) && (dx <= 0))
      {
         // Get points on right of xvalue
         for (int j = start; j < end; j++)
         if (i != j)
         {
            int dx = pts.point[j].x - xvalue;
            if ((0 <= dx) && (dx <= dist))
            {
               // Calculate distance if dy is small enough
               int dy = pts.point[j].y - pts.point[i].y;
               if ((-dist < dy) && (dy < dist))
               {
                  int dx = pts.point[j].x - pts.point[i].x;
                  double temp = sqrt(dx*dx + dy*dy);
                  sqrt_count++;
                  if (temp < dist)
                  {
                     dist = temp;
                     p1 = pts.point[i];
                     p2 = pts.point[j];
                  }
               }
            }
         }
      }
   }
}

//-------------------------------------------------------------------
// Recursive algorithm to find closest pair of points.
//-------------------------------------------------------------------
void closest_pair(PointArray &xpts, int xstart, int xend, 
                  PointArray &ypts, int ystart, int yend,
                  Point &p1, Point &p2, double &dist)
{
   // Zero or one point case
   if (xend - xstart < 2)
   {
      dist = MAX_SIZE;
      return;
   }

   // Two point case
   else if (xend - xstart == 2)
   {
      p1 = xpts.point[xstart];
      p2 = xpts.point[xend-1];
      int dx = p1.x - p2.x;
      int dy = p1.y - p2.y;
      dist = sqrt(dx*dx + dy*dy);
      sqrt_count++;
      return;
   }

   // Split data in half 
   PointArray left, right;
   Point l1, l2, r1, r2;
   double ldist, rdist;
   int xmid, xvalue, yvalue;
   tricky_split(xpts, xstart, xend, ypts, ystart, yend, left, right, xmid, xvalue, yvalue);

   // Make recursive calls
   closest_pair(xpts, xstart, xmid, left, 0, left.size, l1, l2, ldist);
   closest_pair(xpts, xmid, xend, right, 0, right.size, r1, r2, rdist);
   
   // Handle simple merge
   if (ldist < rdist)
   {
      p1 = l1;
      p2 = l2;
      dist = ldist;
   }
   else
   {
      p1 = r1;
      p2 = r2;
      dist = rdist;
   }

   // Handle tricky merge
   tricky_merge(xvalue, ypts, ystart, yend, p1, p2, dist);
}

//-------------------------------------------------------------------
// Main program.
//-------------------------------------------------------------------
int main()
{
   // Initialize point arrays
   for (int size = 2; size <= 1000; size++)
   {
      cout << "SIZE " << size << endl;
      PointArray pts;
      pts.random(size);
      PointArray xpts = pts;
      xpts.xsort();
      PointArray ypts = pts;
      ypts.ysort();
   
      // Brute force closest pair
      Point p1, p2;
      double dist, dist2;
      sqrt_count = 0;
      closest_pair(pts, p1, p2, dist); 
      cout << " p1 " << p1.x << "," << p1.y 
           << " p2 " << p2.x << "," << p2.y 
           << " dist " << dist 
           << " sqrt " << sqrt_count << endl;
   
      // Recursive closest pair
      sqrt_count = 0;
      closest_pair(xpts, 0, size, ypts, 0, size, p1, p2, dist2);
      cout << " p1 " << p1.x << "," << p1.y 
           << " p2 " << p2.x << "," << p2.y 
           << " dist " << dist2 
           << " sqrt " << sqrt_count << endl;
      if (dist != dist2)
         cout << "ERROR FOUND " << size << endl;
   }
}